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The stability ol dark solitons generated by supersonic flow of a Bose-Einstein condensate past an 
obstacle is investigated. It is shown that in the reference frame attached to the obstacle a transition 
occurs at some critical value of the flow velocity from absolute instability of dark solitons to their 
convective instability. This leads to the decay of disturbances of solitons at a fixed distance from 
the obstacle and the formation of effectively stable dark solitons. This phenomenon explains the 
surprising stability of the flow picture that has been observed in numerical simulations. 
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1. Introduction. It is well known that plane dark soli- 
tons are unstable with respect to transverse ("snake") 
perturbations. Such an instability was predicted and 
studied in [TJ [2] for the case of shallow solitons described 
by the Kadomtsev-Petviashvili (KP) equation. The in- 
stability of deep dark solitons in a Bose-Einstein con- 
densate (BEC), described by the Gross-Pitaevskii (GP) 
equation was demonstrated in [3]. These analytical pre- 
dictions were confirmed by both numerical simulations 
and experiments with dark solitons in a BEC (see, e.g., 
[I] and references therein) and with optical solitons (see, 
0). 

It is worth noting that such soliton instability is not 
a property of particular equations, but rather a general 
phenomenon of non-dissipative dynamics. Indeed, as fol- 
lows from [pj , the motion of a ID soliton can be described 
by the semi-classical Newton equation m s d 2 X/dt 2 = F s , 
where m s = 2dE s /d(V 2 ) is the "effective mass" of the 
soliton, E s is its energy, V its velocity, and F s is the 
force acting on the soliton. The crucial point is that in 
all relevant cases E s decreases with increasing velocity 
and hence m s < 0. Let us apply now this equation to 
an element of a soliton deformed along the y axis. If the 
wavelength of the perturbation is long enough, the soli- 
ton can be considered as a surface with surface tension 
E s . Then, according to the Laplace formula, the restor- 
ing surface force is equal to F s = E s /R, where R is the 
radius of curvature (see [7], §61). Let the equation of the 
surface be X = Acos(py). In the linear approximation 
R^ 1 w d 2 X/dy 2 — —p 2 X. This results in an unstable 
dispersion relation for small oscillations: 

u = ±i{E s /\m s \) 1/2 p . (1) 



In the GP equation case for which E s oc (c 2 — V 2 ) 3 ^ 2 , 
where c is the speed of sound, Eq. (Til reduces to the 

1/2 — 

equation uj — ±i [(c 2 — ^ 2 )/3] p obtained in [3]. 

Dark solitons in a uniform BEC are described by the 
GP equation which we write here in standard dimension- 
less units with c = 1 : 

i^t = -|AV»+HV s (2) 

where tp(r, t) is the condensate wave function. The solu- 
tion corresponding to a plane dark soliton moving in the 
x direction was found in [5] in the form (k = y/l - V 2 ) 

ip(x + Vt) = [k tanh (k(x + Vt)) - iV] exp(-it) . (3) 

The depth of a dark soliton depends on its velocity V 
which cannot exceed the sound velocity. The dark soli- 
ton is unstable in the whole range of possible parameters. 
Thus, dark solitons formed by means of density or phase 
engineering |H [DJ are unstable with respect to perturba- 
tions depending on the transverse y and z coordinates. 

However, there exists another possibility to generate 
dark solitons in a BEC. As was shown in [ID], dark soli- 
tons can also be generated by a fast enough flow of a 
BEC past an obstacle as is illustrated in Fig. 1 where 
results of a numerical simulation of a two-dimensional 
flow past a disk-shaped obstacle are shown. We see that 
the Mach cone (imaginary lines drawn from the origin at 
angles ±arcsin(l/M) with respect to the x axis, M = 2 
being the supersonic flow velocity at \x\ — > oo) separates 
regions with wave patterns of different nature. Outside 
the Mach cone there is a stationary pattern arising due 
to waves with a Bogoliubov dispersion law radiated by 
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FIG. 1: Wave pattern generated by the flow of a BEC past an 
obstacle located at the origin of the coordinate system. The 
flow velocity is directed from left to right and corresponds to 
the Mach number M = 2. 



the obstacle. It is possible that these waves have been 
observed in recent experiments I12j . Their theory 
has been developed in [12] [13j [14] . Inside the Mach cone 
the nonlinear waves are located. In case of large enough 
obstacles they form dispersive shock waves considered in 
[T5] , but for obstacles with the size about a healing length 
(about unity in our non-dimensional units), as in Fig. 1, 
just one soliton is formed in each symmetrical "shock" . 
At large enough time of evolution and far enough from 
the obstacle, the profile of a dark soliton tends to a sta- 
tionary state described by the solution ^ of the GP 
equation [TU]. This means that dark solitons generated 
by a fast enough flow of a BEC past an obstacle are sta- 
ble in the reference frame connected with the obstacle 
contrary to the above mentioned instability of dark soli- 
tons in the reference frame where the fluid at infinity is 
at rest. We will see below that this seeming contradiction 
is explained by the fact, that the instability of solitons 
in the obstacle reference frame is "convective" only. This 
means that disturbances are convected from the region at 
finite distance behind the obstacle and, hence, the soli- 
tons there become stable. 

It is also known that at small enough velocity of the 
flow the situation is different. Numerical solutions of 
the GP equation demonstrate that a BEC flow past a 
cylindrical obstacle starts to generate vortices at sub- 
sonic velocity M = 0.45 (see [TH1 [TTJ HB]). However, the 
frequency of vortex generation increases with increasing 
flow velocity. Therefore one can expect that at large 
enough velocity vortices are generated so often that the 
mean distance between them becomes less than their ra- 
dial size so that their separation from each other takes a 
long time which results in the formation of dark oblique 
solitons at a finite distance from the obstacle. 

2. Convective and absolute instabilities. The presence 



of a positive imaginary part in the dispersion law lj — 
Lo{p) of the soliton oscillations does not mean in itself that 
an arbitrary perturbation in a given point will grow with 
time. The stability of the soliton is determined by the 
asymptotic behavior of wave packets built from harmonic 
waves. Let the Fourier transform of the initial density 
perturbation be Sn p . Then the time dependence of the 
perturbation is defined by the equation 

/oo 
5n p e^ py '^ p ^dp . (4) 
-oo 

Hence, there are two possibilities. In one case, the per- 
turbation increases without limit at any fixed value of y. 
This situation corresponds to absolute instability. In the 
other case the packet is carried away by the flow along the 
soliton plane so fast that the perturbation tends to zero 
as t — ► oo at any fixed value of y. This corresponds to 
convective instability. We emphasize that from a practi- 
cal point of view a convective instability does not prevent 
the observation of a finite size soliton in experiments and 
does not violate the stability of numerical simulations. 

The asymptotic behavior of the perturbation Q at 
t — ► oo can be investigated by deformation of the con- 
tour of integration in the complex p plane. The crite- 
rion distinguishing absolute and convective instabilities 
can be formulated as follows (see [19], Ch.VI): (i) in the 
complex p plane the function p(ui) has values lying in 
upper and lower half-planes for Imw = u" — > +oo; (m) 
when we decrease u" then for some values of lu the two 
values of p moving from opposite half-planes coincide in 
the branching point p br of the function p(u>); (Hi) if these 
branching points correspond to values of u} br lying in the 
upper u> complex half-plane, then we have absolute in- 
stability; otherwise the instability is convective. The 
asymptotic behavior of the perturbation at t — > oo is 
Sn(t, y) oc e^P^-^t) w h ere p br = p(uj br ). 

This investigation demands knowledge of the disper- 
sion law uj(p) for all values p so that knowledge of just the 
long wave-length approximation ([!]) is not sufficient for 
this aim. Before going to concrete calculation, we have 
to make a following general remark. The solution (|3| de- 
scribes a soliton in the absence of a flow along its plane. 
In our problem the soliton is created by the flow moving 
with respect the obstacle with velocity M{= Mc). Let 
8 be the angle between a vector normal to the soliton 
plane and the direction of the flow. Then the normal 
component of velocity is V = M cos and the flow has 
a component u — Msm6 along the soliton plane. (A 
solution of Q corresponding to this "oblique" soliton 
has been obtained in [TU].) Since we want to investigate 
the stability of the soliton in the "laboratory" coordinate 
frame with the obstacle at rest, we have to transform the 
dispersion law luo(p) of the oscillation of a "direct" soli- 
ton Q to the "laboratory" reference frame which means 
that ujq(p) should be replaced by w(p) = up + luo(p)- 
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3. Shallow dark solitons. The investigation of the na- 
ture of instability of the soliton ([3| is a quite complicated 
problem. As a first step, we will study a soliton with a 
small amplitude, which moves with the near-sonic veloc- 
ity V, 1 — V -C 1. In this case the dispersion law can be 
found analytically. As is known, in the small amplitude 
limit the GP equation can be reduced to the KP equation 
(see, e.g. [3]) 



-,nn x 



2 vv 



(5) 



where n = 1 + n and we used the same notations t, x, y 
for the new variables. 

This equation has the soliton solution 



n(x,t) = — s/ cosh 2 [y / i(a; + Vt)] 



(6) 



where s — 2(1 — V). Since in the small amplitude limit we 
have 1 — V 2 = 2(1 — V), this solution is an approximation 
to the exact solution ((3| . Small oscillations of the soliton 
^ were studied in [2l~l20] where, for waves propagating 
with wave number p along soliton, the spectrum Wo(p) — 
iT(jp) was obtained, which after transformation to the 
"obstacle" frame takes the form 



uj = ui(p) = up+iT(p) = up+i(p/ V3)y s — 2p/V3 • (7) 

The above criterion can be easily applied to the dis- 
persion relation Indeed, at uj" — > +oo this relation 
yields three values of p one of which (p tx — (w") 3 / 2 ) lies 
on the real axis, and the other two (p oc (ci/') 3 / 2 e ±7 ™/ 3 ) 
are located in the upper and lower p half-planes. Then, in 
the branching point we have du(p) /dp — and this equa- 
tion gives after a simple calculation the critical value of 
p: 



Pb, 



= [s + u 2 - u\/u 2 - s]/V3. 



(8) 



Substitution of this value into Eq. ^ yields the critical 
value of uj: 



^br — Pbr(V U 2 — S + 2ll)/3. 



(9) 



For u 2 < s we get uj% r = ((s - m 2 )/3) 3/2 > which 
means the instability is absolute. However, pi, r and u>b r 
have real values for u 2 > s and in this case the insta- 
bility is convective. Thus, the transition from absolute 
to convective instability takes place at u 2 > s, that is at 
M 2 > 1. The shallow oblique solitons are convectively 
unstable therefore for supersonic flow of a BEC past an 
obstacle. 

4- Convective instability criterion for deep dark soli- 
tons. In the case of deep dark solitons there is no analyt- 
ical expression for the full instability spectrum, although 
its asymptotic expressions are known: for small p -C q 
we have (see [5]) 

T 2 (p) = k 2 p 2 /3-(2-k 2 )p 3 /(3V3), (10) 




FIG. 2: Growth rate approximated by Eq. ( 13 1 with f(p) 
given by the interpolation formula as a function of wave vec- 
tor p of a harmonic wave propagating along the oblique dark 
solitons with V = 0.7, V = 0.5, V = 0. The exact dependence 
found numerically for V = is shown by dots for illustration. 



where k 2 = 1 — V 2 and 



-(1 + V 2 ) + -1/2 + 1 



1/2 



(11) 



Notice that in the limit p — > Eq. (10 1 coincides with 
the general result obtained for the GP equation. For 
q — p <C q we have (see [21] ) 



T 2 (p) = q(q-p)/(3(k), q-p<^q, 



(12) 



where /3(fc) is evaluated numerically. We approximate 
this behavior by an interpolating function 

T 2 ( P ) = f(p)(q-p) (13) 
with f(p) = a(k)p 2 + b(k)p 3 + c{k)p A , where 

a(k) = k 2 /(3q),b(k) = [k 2 /q - (2 - k 2 )/V3]/(3q), 
c(k) = [3//3(fe) - 2fc 2 + q(2 - k 2 )/VS]/(3q 3 ), 

and P(k) can also be approximated as (3{k) = 3/[(l + 
ak 2 )k 2 ], a = 0.596. This interpolating expression 
yields the dependence of T on p shown in Fig. 2 for several 
values of V. 

On the contrary to the shallow soliton case, an expres- 
sion for the critical flow velocity M cr cannot be given in 
explicit form. Therefore we shall follow Sturrock's [22 
formulation of the criterion and notice that the func- 
tion p = p{uj) determined implicitly by the equation 
uj(p) = up + iT(jp,V), changes its behavior at 
for u < u cr it is represented in the complex p plane by 
disconnected curves, whereas for u > u cr these curves 
are connected with each other. Therefore a "space- like" 
Fourier representation of a wave packet of disturbance 
cannot be transformed to its "time-like" representation 
for u < u cr and in this case the instability is absolute. 
For u > u cr we can deform the contour p so that a 
single- valued function uj(p) is defined for all p > and 
hence a "space-like" packet is also a "time-like" packet 
which means that the instability is convective (see [55] 
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FIG. 3: The boundary between regions of convective and ab- 
solute instabilities in the plane of parameters (tan#, M cr ). 
Below this curve the instability is absolute and above it it is 
convective. Dots correspond to the transition points calcu- 
lated from a numerical solution of the eigenvalue problem for 
the instability spectrum. 



for details) . At the branching point of p = p(oj) we have 
dw/dp = 0, that is the values pb r satisfy the equation 



= -idT/dp= -iT'(p,V), 



(14) 



where T(p, V) has either real or purely imaginary values 
for real p. Hence, the critical value of u corresponds 
to the appearance of a double zero p\, r of the equation 



( 14 1 on the real axis of p as it takes place in the above 
considered case of shallow solitons, that is dpt, r /du = oo 
at (In other words, Pb r has here a branching 



point as a function of u.) Differentiation of Eq. ( 14 1 with 
respect to u then gives the equation d 2 T / dp 2 \ ri _ 
for the corresponding critical value p cr (V) 
of p cr , found in this way, into squared equation ( 14 1 yields 
with account of |7| the equation 



= 

P=Pcr 

Substitution 



Ml sin 2 



e=-[T' (p cr (M cr cos 8),M cr cos 8)} (15) 



which determines implicitly the critical value M cr of the 
transition from absolute instability to a convective one as 
a function of the angle 8. For T(p) given by Eq. Q these 
formulae reproduce the above results derived for shallow 
solitons. The plot of M cr (tan6*) calculated for the insta- 
bility spectrum defined by the interpolation formula ( 13 ) 
is shown in Fig. 3 by a solid line. As we see, the curve 
obtained with the use of the interpolation formula agrees 
well enough with the results derived from the numerically 
calculated instability spectrum. 

In experiment [TT] , the obstacle was created by a laser 
beam. Hence the flow of BEC past a cylindrical obsta- 
cle was actually studied. Our results can be applied to 
this case without any changes due to the isotropy of the 
growth rate spectrum T(p), p = (p 2 + p 2 ) 1 ^ 2 , in mo- 
mentum space. Due to this isotropy, the unstable modes 
grow up with the same increments in all directions and 
if the addition of convection can remove a wave packet 
built from unstable modes from a region finite in the x 
direction, then it removes it from a region finite in the y 
direction, too. 



5. Conclusion. Thus, for M < 1 oblique solitons are 
absolutely unstable with respect to the transverse insta- 
bility leading to their decay to vortex pairs of opposite 
polarities. However, for supersonic flow M > 1 a re- 
gion of angles 8 appears where oblique solitons are only 
convectively unstable which provides the possibility of 
formation of oblique solitons generated by the flow of a 
BEC past an obstacle. For M greater than some maximal 
value M max = 1.46, oblique solitons are only convectively 
unstable for any angle 8 inside the Mach cone. These re- 
sults explain the stability of oblique solitons observed in 
numerical simulations [TU] ■ We suppose that this method 
of generation of dark solitons in a BEC presents new pos- 
sibilities for their experimental study. 
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